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ABSTRACT 

Bearing in mind the application to core-collapse supernovae, we study nonlinear properties of the 
magneto-rotational instability (MRI) by means of three-dimensional simulations in the framework 
of a local shearing box approximation. By changing systematically the shear rates that symbolize 
the degree of differential rotation in nascent proto- neutron stars (PNSs), we derive a scaling relation 
between the turbulent stress sustained by the MRI and the shear-vorticity ratio. Our parametric 
survey shows a power- law scaling between the turbulent stress (((u>tot))) and the shear-vorticity ratio 



{g q 



as 



{{wtot)) oc g s q with its index 5 



0.5. The MRI-amphhed magnetic energy has a similar 



scaling relative to the turbulent stress, while the Maxwell stress has slightly smaller power-law index 
(~ 0.36). By modeling the effect of viscous heating rates due to the MRI turbulence, we show that 
the stronger magnetic fields or the larger shear rates initially imposed lead to the higher dissipation 
rates. For a rapidly rotating PNS with the spin period in milliseconds and with strong magnetic fields 
of 10 15 G, the energy dissipation rate is estimated to exceed 10 51 erg sec -1 . Our results suggest that 
the conventional magnetohydrodynamic (MHD) mechanism of core-collapse supernovae is likely to be 
affected by the MRI-driven turbulence, which we speculate, on one hand, could harm the MHD-driven 
explosions due to the dissipation of the shear rotational energy at the PNS surface, on the other hand 
the energy deposition there might be potentially favorable for the working of the neutrino-heating 
mechanism. 

Subject headings: instabilities — rotation — supernova: magnetic fields — stars: proto-neutron stars 



1. INTRODUCTION 

Numerical simulations of magnetohydrodynamic 
(MHD) stellar explosions have started already in the 
early 1970's shortly af t er the discovery of pulsars 
dLeBlanc fc Wilson! 119701: iBisnovatvi-Kqgan et all 1 19761 : 
iMiiller fc HillebrandtJ 119791: ISvmbalistvl ll984ft How- 
ever, it is rather only recently that the MHD studies 
come back to the front-end topics in the supernova 

research followe d by a numb er of extensive MHD 

simul ati ons (e.g.. lArdelian et all 120001: lYamada fc Sawai) 
20041: IKotake et al l l2004all9; lObergaulinger et all 
2006allbT iBurrows et all |2007t iCerda-Duran et al 



2007t iTakiwaki et all 120091: IScheidegger et all 



2010; 



2011 



Takiwaki fc Kotakd 20111 : lObergaulinger fc Jankal . 
and IKotake et al" I 12006T l2012allbD for recent reviews 
Main reasons for this activity are obse rvations in- 
dicat in g ve ry asymmetric explosions (Wa ng~et al.l 
120011 I20021). and the interpretation of magnetars 



(jDuncan fc Thompson! 1 1992t lLattimer fc Prakashi [20071 



and gamma-ray bursts (e.g., IWooslev fc Heger] 1200 
lYoon fc Langerl l2005h as a possible outcome of the 
magnetorotational core-collapse of massive stars. 

The MHD mechanism of stellar explosions relies on the 
extraction of rotational free energy of collapsing progen- 
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itor core via magnetic fields. Hence a high angular mo- 
mentum of t he core is precon ditioned for facilitating the 
mechanism (|Meier et al.lll976l ). Given (a rapid) rotation 
of the pre-collapse core, there are at least two ways to 
amplify the initial magnetic fields to a dynamically im- 
portant strength, namely by the field wrapping by means 
of differential rotation that naturally develops in the col- 
lapsing core, and b y the so-called magneto- rotational in- 
st ability (MRI, se e lgalbus fc Hawlevlll998l ) . 

lAkivama et al.l (|2003l ) were the first to point out that 
the interfaces surrounding the nascent protoneutron stars 
(PNSs) quite generally satisfy the instability criteria for 
the MRI. Therefore any seed magnetic fields can be am- 
plified exponentially in the differentially rotating layers, 
much faster than the linear amplification due to the field 
wrapping. After the MRI enters to the saturated state, 
the field strength might reach ~ 10 15 ~ 16 G, which is high 
enough to affect the supernova dynamics. Not only to 
amplify the magnetic fields, the MRI plays a crucial role 
also in operating the MHD turbulence (seelHawley et al.l 
fl995t iBaibus fc Hawlevl[l998t IMasada et al.ll2006T) . The 
turbulent viscosity sustained by the MRI can convert a 
fraction of the shear r otational energy to the t hermal en- 
ergy of the system. iThompson et al.l (|2005l ) suggested 
that the additional energy input by turbulent viscous 
heating can help the neutrino-driven supernova explo- 
sion. Followed by the exponential field amplification and 
the additional heating, a natural outcome of the magne- 
torotational core-collapse may be the formation of ener- 
getic bipolar explosions, which might be obser ved as the 
so-called hypernovae (see iTanaka et al.l 120091 and refer- 
ences therein). 

Here it is noted that bipolar explosions obtained in 
the previous MHD supernova simulations mentioned 
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above are not driven by the MRI, but predominantly 
by the field wrapping, assuming very strong pre-collapse 
magnetic fields (> 10 12 G) in general. The growth 
rate of MRI-unstable modes depends on the product 
of the initial field strength and the wavenumber of 
the mode. In the case of the canonical initial fields 
(~ 10 9 G) as predicted by recent stellar evolution mod- 
els (|Heger et al.l I2005T) , the fastest growing modes are 
estimated to be at most a few meters in the collaps- 
ing iron core (|Ober gaulingcr ct al. 2009). Unfortunately 
however, it is still computationally too expensive to re- 
solve those scales in the global MHD simulations, typi- 
cally more than two orders-of-magnitudes smaller than 
the typical finest grid size. To reveal the nature of the 
MRI, local simulations focusing on a small part of the 
MRI-unstable regions are expected to be quite useful as 
traditionally studied in the context of accretion disks (see 
iBalbus fc Haw!evl[l998l) 

lObergaulinger et al.1 (|2009l) were the first to report 
their numerical simulations to study the linear growth 
and nonlinear properties of the MRI in the supernova 
environment. To ease a drawback of the local shearing 
box simulation, they employed the shearing disk bound- 
ary conditions by which the global radial density in the 
vicinity of the equatorial region in the supernova core 
can be taken into account. By performing such a semi- 
global simulation systematically in two (2D) and three 
dimensions (3D), they derived scaling laws for the ter- 
m ination of the linear g rowth of the MRI. As estimated 
in lAkivama et a l. (2003), the MRI was shown to amplify 
the seed fields exceeding 10 15 G. These important findings 
may open several questions that motivate us to join in 
this effort, such as how the non-linear properties as well 
as the scaling laws could be changed by other parameters 
yet unexplored in the supernova context (that we are go- 
ing to explain in the next paragraph), and whether the 
viscous heating maintained by the MRI turbulence could 
or could notaffect the supernova mechanism. 

A fundamental and long-lasting issue regarding the 
MRI itself is to understand the features of the MRI in 
the non-linear phase and to specify which physical quan- 
tities determine the saturation levels of the MRI. So far, 
extensive efforts have been paid to measure parameter 
dependence of the MRI-sustained turbulence by means 
of numerical simulations. It was reported that the am- 
plitude of the turbulent stress maintained by the MRI 
depends on the net value of the initial magnetic field 
and the gas pressure of the system, that is 

w tot oc £?«P C , (1) 

where £ and C are the power law indices (see also Black- 
man et al. 2007 for the scaling relation between the tur- 
bulent stress and the plasma beta). Hawley et al. (1995) 
found that the saturation amplitude depends on the field 
strength when the system is penetrated initially by a uni- 
form vertical field (w to t °c B), but is independent of the 
initial field strength if there is no net magnetic flux in 
the system (see also Sano et al. 2004). However, in the 
density stratified system, it is recently suggested that 
the turbulent stress increases almost linearly with mag- 
netic energy of the net vertical field (tUtot fx B 2 ; Suzuki 
et al. 2010; Okuzumi & Hirose 2011). A weak depen- 
dence of the gas pressure on the MRI turbulence was 
found by Sano et al. (2004): ( ~ 1/4 for models with 



initial net-zero magnetic flux and £ ~ 1/6 in the sys- 
tem penetrated by uniform magnetic field. The physical 
mechanisms responsible for these parameterdependences 
remain an issue under considerable debate. 

In the stellar interior condition, such as in the super- 
nova core, it should be important to understand how 
the degree of differential rotation of the PNS could af- 
fect the nonlinear properties of the MRI. This is because 
the shear rate at the PNS surface (q = — d\nil/d\nr 
with SI and r representing the angular velocity and ra- 
dius, respectively) is time- and location- dependent un- 
like accretion disks in which a force balance between the 
centrifugal force and gravity is generally maintained. 

In this work, we study the MRI-driven turbulence 
and their nonlinear properties by performing three- 
dimensional simulations in the framework of a local 
shearing box approximation. By changing systemati- 
cally the shear rates that symbolize the degree of dif- 
ferential rotation in nascent PNSs, we specifically study 
how the nonlinear properties of the MRI-driven turbu- 
lence respond to it. Based on our numerical results, we 
estimate the energy deposition rate due to the viscous 
heating driven by the MRI and discuss its potential im- 
pacts on the supernova mechanism. 

This paper opens with descriptions of the numerical 
methods and the initial settings in § 2. The numerical 
results are presented in § 3. Based on the numerical 
results, we then move on to address possible impacts of 
the MRI on the supernova mechanism in § 4 followed by 
discussions in § 5. We summarize our new findings in § 6. 

2. NUMERICAL METHODS AND INITIAL SETTINGS 

To study nonlinear properties of MRI, MHD equations 
are solved with a MPI-paralleliz ed finite-differenc e code 
that was originally developed by Sa no et all JT998). The 
hydrodynamic part is based o n the second-order Go- 
dunov scheme ()van Leerl Il979f) , which consists of La- 
grangian and remap steps. The exact Riemann solver 
is modified to account for the effect of tangential mag- 
netic fields. The field evolution is calculated with Con- 
sistent MoC-CT method and thus can avoid the numer- 
ical explosive instability appeared in the strong shea r 
layer when the MoC method is adopted ([Clarke! 11996). 
The energy equation is solved in the conservative form. 
The advantages of our scheme are its robustness for 
strong shocks and the satisfact ion of the divergence- 
free constraint of magn etic fields (|Evans fc HawlevHl98l 
IStone fc Normanlll992D . 

We perform a series of 3D compressible MHD simula- 
tion by ado pting the local sheari ng box model described 
in detail by lHawlev et al.l (|T995) • In the shearing box 
model, MHD equations are written in a local Cartesian 
frame of reference (x, y, z) co-rotating with the local 
portion of the stellar interior rotating at the angular ve- 
locity f2 corresponding to a fiducial radius in the cylin- 
drical coordinate R. Then the coordinates are presented 
as x — r — R, y = R(f> — fit, and z. The fundamen- 
tal equations are written in terms of these coordinates 
within a small region surrounding the fiducial radius, in 
Ar « i?, 

^ + v-(H=o, (2) 
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2flxv- 2qn 2 xx , (3) 

(4) 



dB 

~8t 



: V x (v x B) , 



(5) 



where e is the specific internal energy, and the other pa- 
rameters have their usual meanings. The term — 2qil 2 x 
in the momentum equation is the tidal expansion of the 
effective potential with a shear rate q = —dlnQ/dhir. 
Assuming the ideal gas, the pressure is given by P = 
(7 — l)pe. The constant ratio of specific heats 7 is consid- 
ered for simplicity. We choose adiabatic gas with 7 = 5/3 
in this paper. 

Since we focus on the local properties of the instabil- 
ity, we employ, as is described schematically in Figure 1, 
numerical grid representing a small portion of the convec- 
tively stable upper PNS below the neutrinosphere where 
the strong shear is naturally developed when the corc- 
collapsc (Akiyama et al. 2003; Thompson et al. 2005). 
Adopting the differentially rotating matter as an unper- 
turbed state, the azimuthal velocity is given by v y = q£lx 
in the frame co-rotating with the velocity RQ. The radial 
force balance at the initial state is thus achieved between 
the Coriolis force and the tidal force which is the resid- 
ual of the gravitational force mainly balancing with the 
background pressure gradient force. We discuss an appli- 
cability of our numerical results to the high-temperature 
and highly stratified plasma realized in the supernova 
environment and the effects of neutrino viscosity in §5. 

We choose normalization with po — 1, = 10~ 3 , and 
the computational domain has a radial size L x = 4, an 
azimuthal size L y = 4, and a vertical size L z = 1. All 
of the runs use a uniform grid of 256 x 256 x 64 zones. 
The initial field geometry is a uniform vertical magnetic 
field B = P>Qe z (net non-zero flux). We assume, for all 
the models, that the initial gas pressure is Po = 5 x 10~ 7 
and the initial ratio of the gas and magnetic pressures, 
that is the plasma beta is /3q = 3200. The pressure scale 
height is then H p = 1, and is the same as the vertical box 
size. With this normalization, the initial field strength 
is B = 6.26 x 10~ 5 , yielding v A0 = 1.77 x 10~ 5 , where 
vaq = B$/ (Airpo) 1 / 2 is the Alfven speed. Note that these 
non-dimensional parameters can be translated to the di- 
mensional form suitable for supernova environments as 
p = 10 12 gem -3 , L = 10 2 cm, Q, = 10 2 rad sec -1 , 
P = 5 x 10 25 dyn cm -2 , and B = 2.2 x 10 12 G, which 
are summarized in Figure 1. 

The initial characteristic wavelength of the MRI is 
given by Amri = ^V v Azo/^ f° r ideal MHD case 
(|Balbus k Hawlevl 119981) . where 77 = 1/V 4 <7 - <7 2 - Its 
ratio to the vertical box size is thus varied as \mri/L z = 
0.22r) in our models. It is empirically known that, to gain 
a saturation level correctly, the MRI wavelength must be 
resolved at the saturated state by at least six grid zones 
(see Sano et al. 2004) . Here we choose the initial setting 
which satisfies this empirical rule for all the models. 

3. NUMERICAL RESULTS 
3.1. A Fiducial Run 



As a fiducial run, we examine temporal evolutions of 
the MRI in a model with a shear rate q = 0.5. Initial 
perturbations are introduced by giving a random velocity 
perturbations which are taken to have a zero mean value 
with a maximum amplitude of \Sv\/ \J '7P0/ po = 5 x 10~ 3 . 

To characterize the properties of the MRI-driven tur- 
bulence, we pay particular attention to the magnetic en- 
ergy, the Maxwell and Reynolds stresses, and total tur- 
bulent stress defined respectively by, 



E, 



\B\ 



mag 

W M = -BrB^/AlT, 

w R =pv r Sv^, 

Wtot = W M + WR , 



(6) 
(7) 
(8) 
(9) 

where dv^ is the perturbed azimuthal velocity. Note that 
a volume average of these quantities is represented by 
single brackets (like {E mag }) and a time- and volume- 
average is by double brackets (like ((E mag ))). 

The thick, dash-dotted and dashed curves in Figure 2 
show the temporal evolution of volume-averaged mag- 
netic energy, Maxwell and Reynolds stresses normalized 
by the initial gas pressure Po- Note that the horizontal 
axis is normalized by the rotation period t TO t = 2ir/fl = 
6.28 x 10 3 (62.8 [ms] in dimensional form). In order for 
the demonstration of the typical evolution of the MRI, 
four time snapshots of the 3D structure of the radial mag- 
netic field are visualized by means of volume rendering 
method in Figure 3. The color bar indicates the ampli- 
tude of the radial magnetic field, yellow for the positive 
values, and blue for the negative value. Panels (a)-(d) 
correspond to the snapshots at the times t — 6i ro t, 8£ ro t, 
9i ro t, and llt ro t respectively. 

As was described in the previous MRI studies for accre- 
tion disks with Keplerian rotation of q = 1.5, there are 
three typical evolutionary stages observed in our fidu- 
cial run, those are (i) linear exponential growth stage, 
ii) transition stage, and (iii) nonlinear tur bulent stage 
Hawlev k Balbus!ll992tlHawlev et al]|1995l ). Each stage 
is denoted by white, dark gray and light gray shaded re- 
gions in Figure 2. 

In the stage (i), the channel structure of the magnetic 
field exponentially evolves and inversely cascades to the 
larger spatial scale with small structures merging as the 
magnetic field is amplified. This is because the chan- 
nel mode of the MRI is a n exact solution even for the 
nonlinear MHD equation (jGoodman k Xulll994l) . The 
temporal evolution of channel structures for the radial 
magnetic field is shown in the panel (a) of Figure 3. 

Then in the stage (ii), the channel structure of the 
magnetic field is disrupted via the parasitic instability 
and/or magnetic reconnection at the transition stage 
(|Goodman k Xvj|1994t lObergaulinger et al.ll2009D . The 
channel disruption induces a drastic phase-shift from the 
coherent structure to the turbulent tangled structure of 
the field as is illustrated by the panels (b) and (c) of 
Figure 3. The magnetic energy stored in the amplified 
magnetic field is then converted to the thermal energy of 
the system. 

Finally in the stage (iii), the nonlinear turbulent state 
emerges after the channel disruption is maintained long 
enough by the MRI-driven turbulence as is demonstrated 
in the panel (d) of Figure 3. Figure 2 shows that, at 
this stage, the turbulent Maxwell stress dominates the 
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Reynolds stress. It then converts the free energy stored 
in the differential rotation into the thermal energy of the 
system. The time- and volume-averaged total turbulent 
stress is saturated, in the fiducial model, at the level 



((w tot )} ~ 0.03P 



(10) 



where the time average is taken over the saturated state 
from 50t rot to 150t rot . It can be stated that the variation 
of the shear rate does not change the nonlinear properties 
of the MRI qualitatively. 

In the shearing box model adopted in our study, 
such free energy is continuously injected from the radial 
boundaries. This is the reason why the turbulent stage 
powered by the MRI is maintained in our simulation. We 
discuss in § 5 the applicability of our numerical model in 
comparison with th e shea ring disk model employed in 
lObergaulinger et all (|2009t ). 

3.2. Shear Rate Dependence 

In the stellar interior, such as supernova cores, the ro- 
tation profile should be time- and location-dependent un- 
like the typical accretion disk which maintains a quasi- 
Keplerian rotation with the shear rate q ~ 1.5. This 
is because the force balance is mainly achieved between 
the gravity and the pressure gradient force in the stellar 
interior. It should be important to study how the nonlin- 
ear properties of the MRI would respond to the change of 
the shear rate for estimating its impact on the supernova 
dynamics. The shear rates used in our parametric sur- 
vey are listed in Table 1 with a few diagnostic quantities, 
the time- and volume-averaged Maxwell stress, Reynolds 
stress, magnetic energy and the ratio of Maxwell and 
Reynolds stresses, calculated from the numerical data. 

Figure 4 shows the temporal evolution of the turbulent 
stress for the models with different shear rates. The grey, 
red, blue, orange, green and purple curves indicate the 
models with q — 0.2, 0.5, 0.8, 1.2, 1.5 and 1.8 respec- 
tively. Note that the other physical parameters are the 
same as those used in the fiducial run. It is found that 
both the linear growth rate and saturation amplitude of 
the turbulent stress become larger for the model with 
the larger shear rate. This trend can be recognized as 
well from the spatial structures of magnetic fields in the 
saturated state. 

The volume rendered vertical magnetic field is visual- 
ized for the models with (a) q — 0.5 and (b) q = 1.5 
in Figure 5. The time snapshot at t = 100i ro t is chosen 
for both models. The red color corresponds to the posi- 
tive value of the vertical field component and the blue is 
its negative value. The larger structure of the magnetic 
field is found to be generated for the model with the 
larger shear rate. The other components of the magnetic 
field have similar nonlinear properties. 

Figure 6(a) exhibits the time- and volume- averaged 
Maxwell stress ((wm)) (red circles), Reynolds stress 
((wr)) (blue squares), total stress ((w to t)) (green dia- 
monds), and magnetic energy ((E mag )) (yellow triangles) 
as a function of parameter g q , where g q = q/(2 — q) is 
the ratio of vorticity (2 — q)Q, to the shear qQ, (here- 
after called "shear- vorticity ratio", see Abramowicz et 
al. 1996). Note that we take, at the saturated state, the 
temporal average of volume-averaged stresses and mag- 
netic energy from 50i ro t to 150i ro t for models 1,2,3, and 



from 100t r ot to 150i ro t for models 4,5,6,7. The black line, 

1/2 

which is proportional to g q , is plotted as a reference. 

The Maxwell stress dominates the Reynolds stress for 
all the models we surveyed although the fraction of the 
Reynolds stress in the total stress increases with the 
shear- vorticity ratio g q . It should be noted that the 
pure hydrodynamic shear instability develops in the sys- 
tem with the stronger shear of q > 2 (the so-called 
Raylie gh criterion, e . g.. | Cha ndrasckhar 1960). Accord- 
ing to Hawlev et al.l (|1999l ). the hydrodynamic turbu- 
lence is shown to change the properties of MRI turbu- 
lence when the Rayliegh cr iterion is satisfied (see also 
IWorkman fc Armit age 2008). The Reynolds stress then 
exceeds the Maxwell stress. However, since q < 2 
is generally satisfied in the vicinity of the PNSs (e.g . , 
iKotake et all l2004at lObergaulinger et al.l [20Q6al [2009T) . 
the Maxwell stress is expected to dominate over the 
Reynolds stress and plays a major role in the turbulent 
heating process there. 

A power-law relation between the shear-vorticity ratio 
and turbulent stress appears in this figure. The power 
law index 6 of a scaling relation defined by 



((wtot)) oc g s q 



(11) 



is about 1/2 and the best fit is S — 0.44 among all the 
models in Figure 6(a). It would be interesting that the 
magnetic energy has a similar scaling relation to the tur- 

1/2 

bulcnt stress, that is roughly ((E mag )) oc g q (the best 
fit is 8 = 0.46) although the Maxwell stress has a power 
law index of 0.36 (the best fit value) which is a bit smaller 
than that for the total stress. As was well known that the 
maximum growth rate of the MRI is linearly proportional 
to the shear rate q (that is 7 max = qil/2), the scaling law 
of the turbulent stress we found here should be explained 
through the nonlinear properties of the MRI-driven tur- 
bulence. 



3.3. Our Work 



Previous Studies 



Thus far, the dependence of the non-linear proper- 
ties of the MRI on the shear rate has been inves- 
tigated onl y in models whose net magnetic flux is 
set to zero ([Brandenburg et al.lll996l; lAbramowicz et all 
19961: iHawlev et al.l 119991: iZiegler fc Riidigerl 120011 : 
Pessahfc Chanll2008D . These previous studies found the 
stronger dependence of the turbulent stress on the shear- 
vorticity ratio g q (or rather shear rate q itself) which is in 
contrast to the weak dependence (6 — 1/2 see equation 
(fTTj) ) obtained here in the case of the non-zero magnetic 
flux. 

Abramowicz et al. (1996; hereafter ABL96) numer- 
ically evaluated the MRI-sustained turbulent stress for 
different values of the shear-vorticity ratio in the strat- 
ified system with net zero magnetic flux. They found 
that the turbulent stress is roughly represented by a 
linear dependence on the shear-vorticity ratio, that is 
S = 1.0 of equation (|11[) . On the other hand, the nu- 
merical results obtained by Ziegler & Riidiger (2001; 
hereafter ZR01) with a similar numerical setting as that 
adopted in ABL96 suggested that the existence of a rela- 
tion ((wtot)) <x q [i-e., oc 2g q /(g q + 1)] in contrast to the 
finding of ABL96. 

The dependence of the turbulent stress on the shear 
rate was studied comprehensively by Hawley et al. (1999; 
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hereafter HBW99) using unstratified shearing box model 
with net zero magnetic flux . They found that the vor- 
ticity ((2 — g)f2) limits hydrodynamic turbulence and 
strongly reduces the Reynolds stress at around q — 0. In 
contrast, the shear (qfl) promotes turbulence, and thus 
the Maxwell and Reynolds stresses both are enhanced 
with increasing the shear rate. Since the q-dependence 
of the Reynolds stress is much stronger than that of 
the Maxwell stress, the stress ratio diminishes with in- 
creasing q. When extracting the data from Figure 10 
of HBW99, the time- and volume-averaged turbulent 

1 /2 

stress seems to behave as ((u>tot)} oc g q (the best fit 
is 5 = 0.55). 

Liljestrom et al. (2009; hereafter LKKBL09) studied, 
using local shearing box model with net zero magnetic 
flux, the nonlinear properties of MRI-driven turbulence 
with varying shear rate. By fitting numerical data listed 
in Table 1 of LKKBL09, the turbulent stress has strong 
g-dependence like as that of ABL96 and behaves approx- 
imately as ((ui to t)) oc g q (the best fit is 6 = 0.9). 

A theoretical model for the saturation of the Maxwell 
and Reynolds stresses in MRI turbulence was developed 
by Pessah et al. (2006; hereafter PCP06) in light of the 
similarities exhibited by the linear regime of the MRI 
and the turbulent state. On the basis of linear the- 
ory, they formulate a predictor function about the ra- 
tio of the Maxwell and Reynolds stresses at the satu- 
rated state, that is ({wm}}/{{wr}) = (2 + g q )/g q [= 
(A—q)/qm the other form]. LKKBL09 found that the re- 
lation for the stress ratio derived by PCP06 predicts simi- 
lar behavior, but provides a few times smaller ratio than 
that computed from their numerical results. The local 
closure model developed by Ogilvie (2003) rather could 
reproduce numerical results of LKKBL09 quite well. 

For the comparison between our work and the previ- 
ous studies, the ratio of the time- and volume-averaged 
Maxwell and Reynolds stresses is plotted as a func- 
tion of the shear- vorticity ratio g q in Figure 6(b). The 
blue squares correspond to our result, the green crosses 
show the numerical data of HBW99, the orange trian- 
gles are the result of LKKBL09, and the red circles 
present the result of ZR01. A power law relation of 

{{w m)) / {{w r)) oc g q 3 ^ 4 is shown by the solid line for 
the reference. The theoretical model of the stress ratio 
developed by PCP06 is traced by the dashed line. Re- 
member that the previous studies, as was summarized 
above, provided different (^-dependences of the turbu- 
lent stress. The non-zero magnetic flux is imposed only 
in our work though the shearing box approximation is 
adopted in all the models. The stress ratio of HBW99 
is calculated by extracting the data from Figure 10 of 
HBW99. 

The stress ratio decreases with the increase of the 
shear-vorticity ratio g q for all the models. This sug- 
gests that the Reynolds stress has generally a stronger 
g q - dependence than that of the Maxwell stress. As 
LKKBL09 discussed, the model stress ratio of PCP06 
predicts similar behavior, but gives a few times smaller 
value than that obtained by all the numerical models 
especially in the range 0.2 < g q < 10. It is a bit sur- 
prising that the stress ratio is similar among three works 
except HBW99 although they predict different power- 
law relations between the turbulent stress and the shear- 



vorticity ratio. The stress ratio is roughly represented 

by a power-law relation ((wm)) I '({wr}) oc g q 3 ^ 4 in the 
range of 0.5 < q < 1.8. These similar nonlinear charac- 
teristics lied behind the MRI-driven turbulent flows sug- 
gest that the stress ratio might be one of key parameters 
for validating the nonlinear saturation mechanism of the 
MRI. 

4. IMPACTS ON THE SUPERNOVA EXPLOSION 



As taken in iThompson et al.l (|2005h . the a-description 
(|Shakura fc Sunvaevi 119731) conventionally used in the 
context of accretion disks is the simplest way to model 
the turbulent heating. However, it is quite uncertain 
whether the prescription is really applicable to the su- 
pernova problem. Based on our numerical results, we 
hope to reexamine whether the turbulent heating due to 
the MRI could have impacts on the explosion mechanism. 

Combining our result with scaling relations of the tur- 
bulent stress described by equation (1), the heating rate 
maintained by the MRI turbulence can be expressed as, 



= e 



(12) 



Bo) (po) Co) \qoJ vO Q , 

where eo represents the turbulent heating rate at the ref- 
erence state with Bq, Pq, g q o, qo, and f2o. The power 
law indices are positive values depending on the pres- 
ence of the stratification or net magnetic flux, and are 
expected to be £ = 1 ~ 2 (Hawley et al. 1995; Sano et 
al. 2004; Suzuki et al. 2010; Okuzumi & Hirose 2011), 
f = 1/4 ~ 1/6 (Sano et al. 2004), and 5 = 1/2 ~ 1 (sec 
§3.3). 

When choosing the physical parameters adopted in 
the fiducial model as the reference state, that is Bq = 
2.2 x 10 12 G, P = 5 x 10 25 erg cm" 3 , Q = 10 2 rad s"\ 
5q0 = 0-33 and q = 0.5, we can estimate, from equa- 
tion (flu]) , the reference heating rate as e = 0.03P <7o^ — 
10 26 erg cm -3 sec -1 (see §3.1 for the fiducial run). 

We make a crude estimate of the volume in which 
the MRI-driven turbulence becomes active as Vmri = 
4irR 2 h ~ 10 21 cm 3 ^/i 6 , where R 7 = R/10 7 cm is typical 
radius of the PNS normalized by 10 7 cm, and he = h/10 6 
cm is the typical radial thickness of the MRI-active layer 
normalized by 10 6 cm. Then the energy releasing rate 
Lmri = £mriVmri becomes, 

(13) 



7 MRI : 



WR^he erg sec" 



-) ( -) ( ±) - 

Bo) \PoJ \g g oJ \qoJ \V 0/ 
Since we are interested in the MRI in the post-bounce 
stage of core-collapse supernovae, the gas pressure in 
the neutrino opaque upper PNSs is a known parame- 
ter and should be P = 10 31 erg cm -3 . When apply- 
ing the typical gas pressure at the post-bounce phase 
to equation (| 1 3[) with remaining the other parameters, 
the energy releasing rate is slightly enhanced and be- 
comes Lmri — 10 48 erg sec -1 almost independent of (. 
The magnetic field strength and the rotation profile are 
fully unknown parameters in the supernova environment. 
These thus determine whether the MRI-sustained turbu- 
lent heating can assist the supernova explosion or not. 

The energy releasing rate due to the MRI-driven tur- 
bulence Lmri is mapped in Figure 7 as functions of the 
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magnetic field strength B and the shear rate q. Pan- 
els (a) and (b) show the models with the weakest param- 
eter dependence of £ — 1 and S = 1/2, while panels (c) 
and (d) are for the models with the strongest parameter 
dependence of £ = 2 and 5 = 1. The slow rotation of 
ft = fin is assumed in left panels and the fast rotation 
of ft = IOSIq is in right panels. The green, blue and red 
dashed curves correspond to the energy releasing rates 
10 50 , 10 51 , and 10 52 erg sec -1 respectively. The darker 
color gives the higher energy releasing rate. Note that 
the gas pressure is fixed to be 10 31 erg cm -3 and then 
the weak gas pressure dependence of the MRI-driven tur- 
bulence is neglected here for the simplicity. 

The dark shaded parameter space provides the larger 
energy releasing rate than about 10 51 erg sec" 1 which 
should assist the supernova explosion. 5 . The MRI- 
turbulent heating then have a minor effect on the super- 
nova explosion in the unshaded parameter space which 
provides the luminosity Lmri <C 10 51 erg sec" 1 . 

It is fairly obvious that the stronger magnetic field 
or the larger shear rate yields the higher energy re- 
lease due to the MRI, and thus makes a greater con- 
tribution on the supernova explosion. When consider- 
ing the moderate shear of q ~ 1 and angular velocity of 
ft = 0(1O 2 ) rad sec -1 which are plausible for the super- 
nova environment (Ott et al. 2006), a strong magnetic 
field of B > 10 15 G is required for the MRI-assisted su- 
pernova explosion in the case with the weakest parameter 
dependence of £ = 1 and S — 1/2 [Panels (a) and (b)]. 
The magnetic field weaker than 10 13 G can drive the 
MRI-driven turbulence strong enough for assisting the 
explosion when the MRI turbulence depends strongly on 
the physical parameters B and q [see Panels (c) and (d)] . 

The nonlinear MRI studies employing shearing box 
model might suggests that non-canonical post-bounce 
states, such as having strong magnetic field and/or large 
spin rate, must be developed for prompting MRI-assisted 
supernova explosion. It may be tempting to implement 
those phenomenological heating rates as sub-grid model 
into the global MHD supernova simulations to see the 
outcomes. However, before that, it is more important 
to conduct a more extensive nonlinear study of the MRI 
over a wider parameter range and/or using more realis- 
tic numerical setting in order to precisely fix the scaling 
relations. 



5. DISCUSSION 

5.1. Effects of Neutrino Viscosity on MRI 

In the neutrino opaque upper PNSs where we have a 
great interest on in this paper, the neutrino viscosity 
was shown to suppress the growth of the MRI when the 
viscous dissipation timescale of a typical MRI mode be- 
com es shorter than the typical evolution time o f the MRI 
(see IMasada et alj|2007t IMasada fc Sanoll2008D . This is 

5 The typical timescale of neutrino-driven explosions observed 
in the recent supernova simulations is > 300 ms after bounce (e.g . , 
IMarek fc Ja nka 2009; Suwa et al.N20lfl |20"T1 IMiiller et al.H2012l V 
So the energy deposition could be as high as 10 50 erg when the 
MRI-driven turbulence maintains the heating of 10 51 erg sec -1 . 
Since the explosion energy of core-collapse supernovae is esti- 
mated as 10 51 erg, the turbulent heating due to the MRI is 
expected to play an important role in assisting explosion when 
L M m > W 51 erg sec" 1 . 



equivalent to a condition 



R - V A < 1 
fl-MRJ = —77 ^ t, 



(14) 



where -Rmri is the Reynolds number for the MRI, va 
is the Alfven velocity, ft is the angular velocity and v 
is the neutrino viscosity. This can be translated to the 
condition for the magnetic field 



B < B clit = (AirvpVL) 



1/2 



= 3.5 x 10 12 pl^vlfn 



1/2 
2 



[G] 



(15) 



where pi2 is the density normalized by 10 12 g cm" 3 , 
Tu is the temperature normalized by 10 11 K, O2 is the 
angular velocity normalized by 10 2 rad sec -1 , and v is 
the neutrino viscosi ty normalized by 10 10 cm 2 sec -1 (see 
IMasada et al.1120071 for the magnitude of the neutrino vis- 
cosity). Such enormous neutrino viscosity is plausible in 
the region just below the neutrinosphere where the strong 
differential rotation is developed. 

Equation (fT5|) means that if the magnetic field is 
weaker than the critical value even locally, the neutrino 
viscosity suppress the linear growth of the MRI. Lon- 
garetti & Lesur (2010) found that, in the non-resistive 
limit (i.e., highly conducting MHD fluid), the saturation 
level of the MRI does not be affected by the viscosity 
at the nonlinear stage despite the change of the linear 
growth rate. Their results might suggest that the effects 
of neutrino viscosity on the MRI-driven turbulence are 
only secondary in the supernova problem. 

On one hand, when the magnetic field is weaker than 
the critical value B C rit, the linear growth of the MRI 
is suppressed and then the growth time of the fastest 
growing mode of the MRI tMRi becomes 



2(2 -q) 



1/4 



MRI " 



(16) 



(jMasada fc Sand [2008D in contrast to that for the case 
with larger magnetic field than the critical value given 
by £mri = 2/(gft). Figure 8 depicts the growth time of 
the fastest growing mode of the MRI as a function of 
the shear rate for the cases with B > B cr - lt (solid line), 
B = 10" 2 5 crit (dashed line) and B = 10- 4 S crit (dash- 
dotted line). We adopt the spin rate of ft = 100 rad sec -1 
in this figure. 

Considering the spin rate and shear rate plausible for 
the nascent PNSs [ft ~ O(10 2 ) rad sec" 1 and q ~ 1.0], 
the growth time of the MRI becomes longer than 100 
[ms] , which is only marginally comparable to the typical 
timescale of neutrino-driven explo sions observed in some 
recent supe rnova simulations (e.g ., Marck ~fc Jankall2009l : 
ISuwa et al.l 120101: iTakiwaki etaJU l2012f ) when the mag- 
netic field is much weaker than the critical value (dashed 
and dash-dotted curves). This suggests that the weaker 
field can not fully activate the MRI-driven turbulence 
within the dynamical timescale for the supernova explo- 
sion. From the viewpoint of the growth rate of the MRI 
as well as the saturation amplitude, we would require 
the large magnetic field and/or short spin period for the 
MRI-assisted supernova explosion. 

5.2. Validity of Shearing Box Model 

As described in lObergaulinger et al.l ()2009[ ). the shear- 
ing box model can not treat the angular momentum 
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transport adequately. In this respect, the shearing disk 
model is indeed superior to the shearing box model. 
However due to the lack of a simulation technique which 
can appropriately take into account the feed back from 
the global fluid motions, it may not be settled yet (possi- 
bly even by the semi-global simulations) whether or not 
the MRI will really smear out the differential rotation in 
the supernova core. It may be possible that the shear 
energy tapped in the differentially rotating layers can be 
supplied continuously by extracting the thermal or grav- 
itational energies. This is considered to be the case in 
the solar interior. 

It is well known from the helioseismic observation 
that there exists a radial differential rotation in the 
upper part of the convectively stable radiative core 
(the so-called tachocline), where is suitable for the field 
wrapping via fl-e ffect a nd the g rowth of the MRI (see 
iParfrev fe Menoul 120071: iMasadal [20ll . Although the 
field wrapping process and/or MRI might be operated in 
the tachocline region, the radial shear persists secularly 
there. The hydrodynamic angular momentum transport 
processes, such as the meridional circ ulation and the con- 
vective motion (e.g., iRuedigerl [1989) . are considered to 
sustain the differential rotation even if the magnetic ten- 
sion force acts to smear out the differential rotation. In 
the post-bounce supernova core, the convectively stable 
regions are formed between the nascent neutron star and 
the stalled bounce shock, in which th e neutrino c ooling 
dominates over the neutrino heating (|Jankal 120011 ). We 
speculate that the activity of the MRI might be main- 
tained there in the presence of the secular differential 
rotation, which is close to the situation assumed in the 
shearing box simulation. 

To draw a robust conclusion, one should naturally need 
to include the density stratification as well as the effects 
of neutrino heating/cooling which determines the con- 
vective stability. This study is only a prelude to improve 
our modeling according to the long to-do lists one by 
one from now on, to understand the role of MRI on the 
supernova mechanism. 

5.3. Magnetic Field Structure during Saturation 

All our numerical models yield a turbulent, highly- 
tangled, magnetic field structure as the final saturated 
state. The turbulent flow and magnetic field per- 
sist during the saturation, and coherent channel struc- 
tures appear only t ransie ntly. On the contrary, in 
lObergaulinger et "ail ((2009), large-scale coherent fields 
with efficient angular momentum transport emerge after 
the turbulent state and are maintained for some mod- 
els with a uniform initial magnetic field using simulation 
boxes of small radial and azimuthal aspect ratios L x /L z 
and L y /L z . 

The magnetic field structure during saturation would 
be determined by whether flow-driven and current-driven 
parasitic instabilities, which are responsible for the de- 
struc tion of channel solutio ns of the MRI, can grow or 
not. iGoodman fe Xul (|1994|) predicted analytically that 
these parasitic modes require radial and azimuthal wave- 
lengths larger than the vertical wavelengt h of the chan- 
nel solution for being unstable (see also lLatter et al. I 
120091: iPessah fc Goodman! [20091 ). The coherent channel 
structure emerging from developed turbulence is thus ex- 
pected to evolve into the large-scale without the destruc- 



ti on by the parasites in the s mall simulation box. 

lObergaulinger et all ([20091 ) numerically confirmed that 
small radial and azimuthal aspect ratios are required to 
maintain the large-scale structure of coherent magnetic 
fields at the saturated state (see Figs 23 and 24 in their 
paper). A stationary turbulent state with tangled mag- 
netic fields appears even in the shearing disk model when 
the simulation box is large enough. The same trend re- 
garding the magnetic field structure during saturation 
has been reported also in th e shearing box simu lation by 
iBodo et al.l ([2008D (see also iLesaffre et al.ll2009ft . 

As presented in § 2, we chose the computational do- 
main having large aspect ratios of L x /L z — L y /L z = 4 
for all the models. The parasitic instabilities thus can 
evolve without being affected by the box geometry, and 
then destroy coherent channels, yielding the less violent 
turbulent state continuously during t he saturation. It 
should be stressed that, as studied by ISano fc Inutsukal 
([20011 ) . the recurrent formation of large-scale coherent 
fields can be observed in the saturated stage of our shear- 
ing box model when we reduce the simulation domain to 
that with small aspect ratios. 

When the magneto-rotational core-collapse, the MRI- 
active region should be confined in the upper PNS with 
small radial and latitudinal extents of O(10 5 ~10 6 ) cm 
(see § 4). The azimuthal thickness of the MRI- unstable 
region is expected to be significantly larger than the ra- 
dial and latitudinal ones. This geometry would inhibit 
the development of the coherent structure of magnetic 
fields by prompting parasitic instabilities, and lead per- 
sistently to a less violent M RI-driven turbulent state (cf., 
lObergaulinger et"aLll2009l ). The turbulent, highly tan- 
gled magnetic field structure would be suitable for de- 
scribing the MRI- active layer in the supernova cores. 

6. SUMMARY 

We performed a series of 3D compressible MHD sim- 
ulation by adopting the local shearing box model. By 
changing systematically the magnitudes of the shear rate, 
we specifically studied how the nonlinear properties of 
the MRI-driven turbulence are controlled by the shear in 
the system. With applying our numerical results to the 
supernova environment, we examined the impact of the 
MRI on the supernova explosion mechanism. Our main 
findings are summarized as follows. 

1. We could observe, in our fiducial run with q = 0.5, 
three typical evolutionary stages [(i) linear exponential 
growth stage, (ii) transition stage, and (iii) nonlinear tur- 
bulent stage] which are analogous to those found in pre- 
vious works modeling accretion disks. This validates that 
the variation of the shear rate does not change linear and 
nonlinear properties of the MRI qualitatively. The tur- 
bulent stress was saturated at the level ((wtot)) — 0.03Po 
in our fiducial model. 

2. Our parameter survey results in the power-law re- 
lation between the shear-vorticity ratio and turbulent 
stress. The power law index S of ((wtot)) oc g s q is about 
1/2. The MRI- amplified magnetic energy has a simi- 
lar scaling relation to the turbulent stress although the 
Maxwell stress has a power law index of 0.36. 

3. It was found that the stress ratio, defined by 
((wm)}/ ({wr)), decreases with the increase of the shear- 
vorticity ratio. In addition, the stress ratio calculated 
from our numerical results has the similar magnitude and 
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<7q-dependence with those obtained by previous works 
(LKKBL09 and ZR01) despite the different computa- 
tional settings. In the range 0.5 < q < 1.8, the 
stress ratio is roughly fitted by a power law relation 

{{wm))/((w r )) ctg~ 3/4 . 

4. The stronger magnetic field or the larger shear rate 
provides the higher energy release due to the MRI-driven 
turbulence. For a rapidly rotating PNS with the spin 
period in milliseconds and with strong magnetic fields of 
10 15 G, the energy dissipation rate is estimated to exceed 



10 erg sec 



Our results suggest that the conven- 



tional MHD mechanisms of core-collapse supernovae are 
likely to be affected by the MRI-driven turbulence, which 
we speculate, on one hand, could harm the MHD-driven 
explosions due to the dissipation of the shear rotational 
energy at the PNS surface, on the other hand the energy 



deposition there might be potentially favorable for the 
working of the neutrino-heating mechanism. 
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Fig. 1. — Schematic picture of the numerical setting in our calculations. The differentially rotating, neutrino-thick region below the 
neutrinosphere is focused (~ gray shaded region). The temporal evolution of the the local portion of the upper proto- neutron stars is 
studied by the local shearing box model. The physical parameters adopted for our fiducial run are summarized in this figure. 
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Fig. 2. — Temporal evolution of volume-averaged magnetic energy (solid), Maxwell stress (dash-dotted), and Reynolds stress (dashed) 
for the fiducial model with q = 0.5. The vertical and horizontal axes are normalized by initial gas pressure Po and rotation period trot 
respectively. Three typical evolutionary stages, that is (i) linear exponential growth stage, (ii) transition stage, and (iii) nonlinear turbulent 
stage are filled by white, dark gray and light gray shaded colors. 
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Model2 


0.5 


2.86 x 10" 


2 


1.28 x 10" 


-3 


6.57 x 10- 


2 


2.23 x 10 1 


Model3 


0.8 


3.63 x 10" 


2 


3.09 x 10- 


-3 


9.41 x 10- 


2 


1.18 x 10 1 


Model4 


1.0 


3.99 x 10- 


2 


5.15 x 10- 


-3 


1.09 x 10- 


1 


7.78 x 10° 


Model5 


1.2 


5.38 x 10- 


2 


8.49 x 10- 


-3 


1.46 x 10- 


1 


6.33 x 10° 


Model6 


1.5 


6.64 x 10~ 


2 


1.78 x 10- 


-2 


1.94 x 10- 


1 


3.72 x 10° 


Model7 


1.8 


9.63 x 10- 


2 


4.27 x 10- 


-2 


2.96 x 10- 


1 


2.26 x 10° 


TABLE 1 



The shear rates used in our parametric survey are listed in the first column. The second, third, and forth column show 
the time- and volume-averaged maxwell stress ({w m }), reynolds stress ((w r )), and magnetic energy* {(e mag )) normalized 

BY THE INITIAL GAS PRESSURE P = 5 X 10~ 7 . THE RATIO OF MAXWELL AND REYNOLDS STRESSES ((wm)} / {{ w r}) IS IN THE FIFTH 
COLUMN. THE TIME-AVERAGE IS TAKEN OVER THE SATURATED STAGE FROM 50i rot TO 150t rot FOR MODELS 1,2,3, AND FROM 100t rot TO 
150t rot FOR MODELS 4,5,6,7. MODEL 2 CORRESPONDS TO THE FIDUCIAL ONE. 
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Fig. 4. — Temporal evolution of volume-averaged total stress (uitot) normalized by the initial gas pressure Po for the models with different 
shear rates of q = 0.2 (grey), 0.5 (red), 0.8 (blue), 1.2 (orange), 1.5 (green) and 1.8 (purple) respectively. 
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Fig. 5. — 3D volume rendering of the radial magnetic field for the models with (a) q = 0.5, and (b) q = 1.5. The color bar indicates the 
amplitude of the radial magnetic field, that is the red denotes the positive value and the blue is the negative. 




g q = <i/{2-q) g q = q/(2-q) 

Fig. 6. — (a). The volume- and time-averaged Maxwell (red circles), Reynolds (blue squares) and total stresses (green diamonds) as a 
function of the shear- vorticity ratio g q . The yellow triangles represent time- and volume- averaged magnetic energy at the saturated state, 
(b) . The ratio of Maxwell and Reynolds stresses for the models with different initial settings as a function of the shear- vorticity ratio g q . 
The blue squares, green crosses, yellow triangles, and red circles demonstrate the results of our work, Hawley et al. (1999), Liljestrom et 
al. (2009), and Ziegler&Rudiger (2001) respectively. The dashed line traces the prediction from the local model developed by Pesah et al. 

(2006). The solid line gives a power law relation of ( (w m ) ) / ( (wr ) ) oc g q 3//4 . 
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Fig. 7. — The energy releasing rate due to the MRI turbulence Lmri as functions of the magnetic field strength B and the shear rate 
q. Panels (a) and (b) show the models with the weakest parameter dependence of £ = 1 and S = 1/2, while Panels (c) and (d) are for the 
models with the strongest parameter dependence of £ = 2 and 5 = 1. The slow rotation (SI = 100 rad sec -1 ) is assumed in left panels and 
the rapid rotation (f2 = 1000 rad sec -1 ) is in right panels. The green, blue and red dashed curves correspond to the energy releasing rates 
10 50 , 10 51 , and 10 52 erg sec -1 respectively. Note that the darker color gives the higher energy releasing rate. 
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Fig. 8. — The typical linear growth time of the MRI as a function of the shear rate for the cases with B > B C rit (solid line), B = 10 2 B cr i t 



(dashed line) and B = 10 4 -B cr it (dash-dotted line), where B cr i t 
f2 = 100 rad sec -1 in this figure. 
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[G]. We adopt the moderate spin rate of 



